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ABSTRACT 


A matrix partitioning scheme is presented for approximating 
EOF's (empirical orthogonal functions) or eigenvectors of a large 
sample covariance matrix. The data array, a field of measured 
anomalies of some physical variable relative to their time aver- 
ages, is partitioned in either the space domain or the time do- 
main. Eigenvectors and corresponding principal components of the 
smaller dimensioned covariance matrices associated with the par- 
titioned data sets are calculated independently, then joined to 
approximate the eigenstructure of the larger covariance matrix 
associated with the unpartitioned data set. The accuracy of the 
approximation oC (fraction of the total variance in the field) and 
the magnitudes of the largest eigenvalues from the partitioned co- 
variance matrices together determine the number of local EOF's 
and principal components to be joined at any particular level. 

The accuracy and efficiency of the method are demonstrated by 
applying the technique to the space-time distribution of Nimbus-5 
ESMR (Electrically Scanning Microwave Radiometer) sea ice bright- 
ness temperature measurements for a large area of the Antarctic 
(Weddell Sea) and for the time period extending from 30 September 
1973 through 25 May 1975* 

From analysis of the spatial EOF's and their coefficients 
(principal components) maximum and minimum ice extents and the 
times of these extents can be identified. Regions and periods of 
ice growth and decay are identified as well as regions and periods 
of higher changes in growth and decay. The interannual variability 
in the Weddell Polynya between 1973 and 197^ is exhibited by 
the fourth EOF and its principal component. Power spectral analysis 
of the principal components reveal periods which can be related to 
the seasonal cycle of sea ice growth and decay in the Weddell Sea, 
harmonics of this cycle, the cold season (205 days), the warm 
season (l 60 days), and the duration of spring-summer ice removal 
(120 days) as reported in the literature. The first four EOF's and 
their components can be considered the dominant normal modes 
of variation in the ice field, accounting for 85% of the total in- 
formation content in the data (field variance) , 



1.0 INTRODUCTION 


Empirical orthogonal function (EOF) analysis, or the method of 
principal components (also singular value decomposition) is well- 
documented in the literature and has found wide application in many 
areas of applied research, particularly in those areas where there 
are large numbers of measurements. For large data sets (or geophysical 
fields) it becomes increasingly important to find ways of comnressing 
the data while still extracting as much of the information content as 
possible . 

EOF's, which are linear functions of the measurements, provide 
us with an efficient tool in data compression and information extrac- 
tion. The method is efficient in the sense that highly correlated 
fields can be adequately represented by the least number of orthogonal 
functions and corresponding orthogonal coefficients (principal compon- 
ents). The higher the correlation in the measurements, the fewer the 
number of functions and coefficients are required to describe the data 
and explain the variance in the field. This reduces the dimensionality 
of the problem. The latter feature is especially important in the devel 
opment of statistical prediction schemes that rely upon multiple linear 
regression techniques, since the skill and confidence of these schemes 
depend to a large extent upon a priori methods of reducing the number 
of available predictors (Davis, 1976; Barnett and Hasselmann, 1979). 

A property of EOF's which makes them particularly appealing is 
that, unlike conventional orthogonal representations (e.g., the famil- 
iar Fourier decomposition, Tschebycheff polynomials, or spherical har- 
monics), they do not require any predetermined form. Rather, since they 
are derived as the eigenvectors of the covariance matrix between the 
variables (which can be computed for any variable observed on any grid, 
regular— or— no-t-)-,— their— form— depends—di-rectly— on -the— in-terrelationshlps 
existing within the data itself. This feature is especially desirable 
when analyzing fields such as sea surface temperature, pressure, or sea 
ice brightness temperature anomalies, which have no known analytic form 
and depend on complex boundary conditions. Also, the derived eigenvec- 
tors often provide us with a tool for gaining insight into the physical 



interpretation of underlying complex processes within a geophysical 
field. 

The first EOF is that linear combination of the original var- 
iables, which when used as a linear predictor of these variables, ex- 
plains the largest fraction of the total variance. The second, third 
EOF, etc., explain the largest parts of the remaining variance. The 
ordered set of EOF's and their coefficients are frequently referred 
to as the normal (also natural, principal) modes of variation in the 
field. The lower ordered EOF's not only account for the largest 
fraction of the total variance, but their coefficients often show 
large scale frequency variations. The higher ordered EOF's usually 
show small spatial scales with their coefficients being characterized 
by diminished amplitudes and higher frequencies, both of which are 
sometimes associated with noise. 

Early references on principal components can be found in Pearson 
( 1901 ) and Frisch (1929) who were concerned with fitting a line, a 
plane, or in general a subspace to a scatter of points in a higher 
dimensional space. Principal components as applied to random variables 
were introduced by Hotelling (1933) who characterized them by certain 
optimal properties. Since then, they have been characterized by 
slightly different properties (Girschick 1936, Anderson 1958, Kullbach 
1959, Anderson I 963 , Darroch I 965 , Okamoto I 968 ). A good treatment of 
the subject can be found in Kendall (1957), Kshirsagar (1972), Marr- 
iott ( 197 ^), and Kendall (1975)* An excellent paper on the use and in- 
terpretation of principal components in applied research can be found 
in Rao (1964) . This reference contains proofs of a number of optimal 
properties of principal components. One of the more important proper- 
ties states that the eigenvectors of the sample covariance matrix com- 
prise an optimal set of basis functions for a given truncation k. That 
is, for a given truncation k, no other basis set can explain more of 
the average variance (Davis 1976, Appendix A). 

Fukuoka (1951) used principal components in a study of 10-day 
forecasts and Lorenz (1956) was the first investigator to use the term 
"empirical orthogonal functions" (eigenvectors of a sample covariance 



matrix, also characteristic vectors) in connection with statistical 
weather forecasting. Gilman (1957) subsequently applied BOF's to 
30 -day forecasts using data from 40 winters to calculate eigen- 
vectors of temperature over the United States, and pressure over the 
Northern Hemisphere. Sellers (1968) computed eigenvectors for pre- 
cipitation over the western United States for each month of the year 
and Kutzbach (1970) calculated eigenvectors for the January and July 
sea level pressure fields over the Northern Hemisphere. Kidson (1975) 
derived BOF's for temperature, precipitation, and sea level pressure 
in both the Northern and Southern Hemisphere and in the tropics using 
10 years of data. More recently Barnett (1978) studied winter-averaged 
and annually-averaged surface temperature eigenvectors over Northern 
Hemisphere land and ocean areas for the 1950 - 1977 time period. Tren- 
berth (1975) and Davis (1978) used BOF's to study air-sea interactions 
and Walsh and Johnson (1979) analyzed associations between interannual 
atmospheric variability and arctic sea ice extent, 

Bxcellent references for the derivation of BOF's and their 
coefficients as well as procedures for calculating them can be found 
in Kutzbach (I 967 ), Sellers (1968), and Davis (1976). Hotelling (1933)f 
Bartlett (195^) 1 Lawley (1956), Buell (1978), Preisendorf er and. Bar- 
nett ( 1977 * 1978 ) discuss significance tests for BOF's, and Craddock 
and Flood (I 969 ), Craddock and Flintoff (1970), and Rinne and Jarven- 
oja ( 1979 ) outline procedures for truncating the BOF series. 

Buell ( 1972 ) has a very interesting discussion on the integral 
representation of the eigenvalue problem and compares it with the 
matrix formulation. He states "...the integral representation is 
appropriate for meteorological problems since additional considera- 
tions based on the properties of a continuum are possible, necessary, 
and desirable." In two or more dimensions the shape of the boundary 
that defines the area of integration must be taken into account. If 
the area of each grid is the same, the two formulations are equivalent. 
This is essentially true for the data set discussed in the present 
paper. The shape of the boundary is also pursued in Buell (1975» 1979). 

Craddock (1973) discusses the reduction of the dimensionality of 
the problem of long range weather forecasting and mentions significant 
advances in the application of eigenvector techniques in meteorology. 



North, Bell, and Cahalan (I 982 ) give an excellent review of 
EOF's, focusing attention on the necessary weighting factors for 
gridded data and the sampling errors incurred when too small a 
sample is available. A rule of thumb is presented for indicating 
when an EOF is likely to be subject to large sampling errors. 

The number of applications of EOF analysis and the list of 
parameters which have been expanded in terms of EOF's is quite 
extensive. Atmospheric parameters which have been expanded in- 
clude pressure (geopotential), temperature, humidity, water-vapor 
mixing ratio, rainfall, wind velocity components, and ozone con- 
tent (Obukhov i 960 , RukhovetS 19^3 » Holmstrom 19^3 » Grimmer 1963* 
Garrilin 1965» Popov 1965» Koprova and Malkevich 1965* Mateer 1965* 
Marchuk I 965 * Wark and Fleming I 966 , Alishouse et. al 196?, Kutz- 
bach 1967 , Sellers 1968, Yakovleva et. al I 968 , Sellers and Yarger 
1969f Craddock and Flood ± 969 , Craddock and Flintoff 1970, Bodin 
197 ^, Kidson 1975, Weare 1976, 1979, 1982). Spurrell (I 963 ) has 
applied EOF's to metallurgical data, Choi ( 1967 ) to the analysis 
of seismic data. Smith et. al (1972, 1976) to retrieve atmospheric 
parameters from spectral radiance measurements, and Mueller (1976) 
has used EOF/principal component analysis for ocean color spectra, 
Jalickee and Klepczynski (1977) applied singular value decomposit- 
ion (principal component analysis) in the compaction of navigation 
tables (see also Good I 969 and Golub and Reinsch 1970 in connection 
with singular value decomposition) . Steyaert et. al (1978) have 
used EOF's of sea level pressure as predictors of wheat yields in 
North America and the Soviet Union. Barnett and Preisendorf er (1978) 
have formulated eigenvectors of several variables in "key regions", 
identified through a filtering process, in order to analyze climatic 
predictability. And, Walsh and Johnson (1979) have studied Arctic 
sea ice data from the 1953-1977 period, and used EOF's to identify 
the major spatial and temporal scales of ice fluctuations within the 
25 -year period. 

Yakovleva et. al (I 968 ) present a joining method for combining 
spatial EOF's and their corresponding time coefficients from two 
distinct spatial regions in order to approximate the "global" EOF 



structure for the larger domain. This technique is applied in 
Rinne and Karhila (1979) and in this paper. 

The purpose of this paper is two- fold: (1) to present a 
matrix partitioning scheme for approximating EOF's associated 
with large sample covariance matrices at a fixed tolerance 
level (or fixed percentage of the total sample variance), and; 

(2) to apply the method to the EOF/principal component analysis 
of the Nimhus-5 TSMR (Electrically Scanning Microwave Radiometer) 
brightness temperature measurements (1.55 cm) of Antarctic sea 
ice for a large area of the south pole (Weddell Sea region) and 
for the time period from 30 September 1973 through 25 May 1975* 

In addition, a physical interpretation of the first four spatial 
EOF's and their corresponding time coefficients or principal 
components (which can be considered the dominant normal modes of 
variation within the ice field) will be given. 

Section 2.0 is devoted to a brief description of the relevant 
mathematical theory associated with the BOF/principal component 
decomposition. In Section 2.1 the concept of matrix partitioning 
(in either the spatial domain or the time domain) is introduced. In 
Section 3»0» P^'cvious analyses of ESMR data are cited. The relation- 
ships between brightness temperature, emissivity, and sea ice con- 
centration are shown and discussed in Section 3-1* Section 4 con- 
tains a description of the particular data set analyzed, the ESMR 
instrument (spatial and temporal resolution), and the results of an 
EOF decomposition of the field and application of the partitioning 
method to the data set. The conclusions are found in Section 5» 



2.0 MATHEMATICAL THEORY 


In the following V is a real (p x q) data matrix representing 
the scalar field of measured anomalies of some physical variable rel- 
ative to their time averages (row means removed) . The number of 
spatial points is p and the number of time points is q. V has rank 
r ^ q - p. 

The complete singular value decomposition of V (Good 19^9 i Golub 
and Reinsch 1970) or the expansion of V in terms of empirical orthog- 
onal functions (EOF's) and their associated coefficients (principal 
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components) is given by ' ' 
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where the orthonormal vectors e. (columns of E) are the left 
singular vectors of V or spatial EOF's (eigenvectors of VV or of 
the sample covariance matrix^ = VV'/(q-l)), the orthonormal vectors 
f- (columns of F) are the right singular vectors of V or time EOF's 
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(eigenvectors of V'V), the orthogonal vectors y. (row vectors of Y 

J 

normalized to d.) are the principal components in the spatial domain 

J 

corresponding to the e. vectors (functions of time), the orthogonal 

J 

vectors z- (row vectors of Z also normalized to d.) are the principal 
J J 

components in the time domain corresponding to the f^ vectors 

(functions of space) , and the d® are the singular values of V or the 

J 

non-negative square roots of the eigenvalues of VV (or of V'V) 
with 
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(1) x' is the transpose of x. 
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= diag(d| . . . d^) 


E'E = F'F = FF' = I 


YY' = ZZ’ = D = diag(d^...d ) 


Since the rank of V is r “ *•' “ 


Using (1) and the orthonormal property of the e. and f . 

J _ J 

vectors in (2) the principal components in the space domain 
and time domain are, respectively, 




d^e'. 
J J 
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From (2) and (3) » y-;y^ ~ “ d-i ■; « that 

the principal components are orthogonal functions of the data V, 
ordered in such a way that the first component y^ or ( or 
equivalently e^ or f^ ) corresponding to d^ accounts for the 
largest fraction of the trace of VV (or of V'V), the second 
component accounts for the second largest fraction of the 
trace, and the jth component accounts for 


(dj/T^) X 100?^ 


of the trace T where 
r 


= dl + . . . + dq 


From (1) the measurement vector v. is given by 
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The rank of V,V',VV, and V'V is the same (Noble 19^9) 


( 7 ) 



( 5 ) 


q 

k=l 


where the coefficient y, . associated with the eigenvector 

K j 

e^ for the jWi time point is the element of y^^ and repre- 

sents the amplitude at time j of the normalized (unit length) 
spatial EOF pattern e, . 


The sum of squares of the coefficient associated with the 

q ? . 

kth eigenvector, 5^ y, . = d, is equal to (q-1) times the total 
- kj k 


variance explained by the k th eigenvector or EOF. The total 
variance in the field is T^(q-l). 


From (1) we can also write the scalar representation of 
the field as 


q 

''ij = 

k=l 


(6) 


which shows that the contribution to the field by the k th spatial 
EOF at the i th spatial point and time point is given by the 

product of the i^ element of e^^ and the element of the assoc- 

iated time coefficient or principal component yj^. 

Since the elements of y^^ = ( » • • • > ) from (3) are the 

projections of the measured anomaly vectors v- (j = l,2,...,q) ontb 

J 

ej^ (see Appendix A for a geometrical interpretation of principal 
components), the yj^ are time-dependent amplitude functions which 
modulate the spatial EOF patterns e^^ and describe the temporal 
variation of the field about the mean in the direction of e-^. For 
fixed k the e^^, being functions of space only, describe the spatial 
variation of the field about the mean in this direction, i.e., the 
spatial distribution of covariance in the field. The ordered spatial 
EOF's and their associated time coefficients (or principal compon- 
ents) are sometimes referred to as the normal (also natural or 
principal) modes of variation in the field. 


From (1) and (2), V V = Y'Y, from which we obtain 


( 8 ) 



( 7 ) 



and therefore 

(y^j/vjv^) X 100 (8) 


represents the fraction in percent of the total sum of squares 
of the 0 ^ measured anomaly vector (v.) accounted for by the k th 

cl 

EOF e^ (see Sellers I 968 ) or equivalently the percent contribution 
of the k th principal component y^ to the spatially averaged mean 

square anomaly at the time point, v'.v. . 

J J 

In Appendix B it is shown that for fixed k - r, (V-V) has 
minimum norm, where V is the approximation to V using the first k 
terms in (1). Thus, the expansion of V in (1) is optimum in the 
sense of least squares. Each successive term in the expansion more 
closely approximates V with the complete reconstruction of the 
data being accomplished when k = r. For highly correlated data 
arrays or fields the first few EOF's will account for the largest 
fraction of the total information content in the data (trace of 
VV ) . Therefore, representing the field using these few terms 
leads to a reduction in the dimensionality of the problem. It is 
in this sense that the singular value decomposition or EOF/princi- 
pal component expansion of V provides us with the most efficient 
method for data compression (see Jalickee and Klepczynski 1977)* 


It is shown in Appendix B that VV and V V have the same non- 
zero eigenvalues and that if e- is an orthonormal eigenvector of 

c) 

VV corresponding to d . ^ 0 with prindipal component y. , then 

_ i D w 

f. = d.®y'. is an orthonormal eigenvector of VV also corresponding 
3 J J 

to d. + 0. Also, if f. is an orthonormal eigenvector of VV corr- 
J J 

esponding to d . ^ 0 with principal component z., then e. = 

J t) J J u 

Is an orthonormal eigenvector of VV corresponding to d.. The 

lI 

duality between the principal components and eigenvectors in the 
two domains is evident from the relations in (3)* This result can 
afford us savings in both computer storage requirements and compu- 
tation time in the calculation of eigenvalues and eigenvectors, 
especially if one of the dimensions of V is much larger than the 


( 9 ) 



other. The eigenstructure of the smaller dimensioned matrix (either 
W or V'V) can be calculated first. Then, from this structure the 
eigenvectors of the larger dimensioned matrix can be easily obtained. 

If both dimensions of V are large, V can be partitioned into a 
number of smaller subarrays and the eigenstructure of the smaller 
dimensioned matrices - either or - calculated. Then, by 

properly joining eigenvectors and principal components from each 
partition, the eigenstructure of VV (or of V'V) can be obtained to 
any required accuracy o^(0 -o< 1) where represents a fraction of 
the trace T^. The number of vectors and components to be joined at 
any particular level is a function of or', T^, and the magnitudes of 
the largest eigenvalues of V^V|. The method for joining vectors from 
various subdivisions of the data has been presented before by Yak- 
ovleva et. al (1968) and also by Rinne and Karhila (1979). In appli- 
cation of the joining procedure these authors have partitioned the 
data in the spatial domain, joining together a fixed number of 
spatial EOF's from each partition. In the present analysis V is par- 
titioned in both the space and time domains and only that number of 
vectors (and components) are joined at each level which are suffic- 
ient to insure that 100«^^ of the trace T^ is accounted for. 



2.1 MATRIX PARTITIONING AND JOINING VECTORS 


Partitioning of the basic data array V into a number of 
subarrays with the subsequent joining of local EOF's from 
these subdivisions offers two advantages. First, the procedure 
provides us with an efficient algorithm for approximating the 
eigenstructure of VV (or of V'V) at a specified tolerance 
leveled. Second, joining local EOF's from these particular 
subdivisions enables us to gain an insight into the global EOF 
structure through analysis of the joining vectors. The joining 
vectors (eigenvectors of the covariance matrix of principal 
components from the various subdivisions) are essentially 
weighting functions with the elements of each vector weighting 
specific local EOF's to the global EOF associated with that 
particular vector. Since the elements range in value from -1 
to +1 joining vectors can also be viewed as correlation funct- 
ions. In this sense each element gives us a measure of the 
correlation between the local EOF and the global EOF. These in- 
terpretations of the joining vectors are useful in interpreting 
global EOF's. 

Sets of EOF's associated with the various subdivisions of 
the data may be joined together in a number of ways. The pro- 
cedure for joining M sets in one operation, is an extension of 
the procedure for joining two sets in one operation ( this 
latter procedure is described in Yakovleva et. al 19^8 and in 
Rinne and Karhila 1979) and is shown in Appendix C. Appendix D 
contains a general procedure for joining where there are L 
partition levels, matrices and groups of matrices at the 
Ith level (1 = 1 ,2, . . . ,L) . 

For purposes of this discussion we will consider partit- 
ioning in the spatial domain only. To apply the procedure of 
joining in the time domain the transpose of V, or V , can be 
used. The algebraic equations and manipulations which apply 
for V will then apply for V . Thus, let us partition V in (1) 
into M submatrices each of dimension (p^ x q) and rank 

r^ ^ min(p^,q) with + • • • + Pjj = P and r^ + . . . + r^^ - r. 
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(9) 


From Appendix C, V can be written as 


V = EY = E(AD®H' ) 

= GW = + ... + gj^w^ (10) 


where 





The columns of E^ in E are the local spatial BOP* s for the 
i th subdivision and the rows of in Y are the corresponding 
principal components. Since the rows of Y are correlatedjY is 
expressed in terms of its singular value decomposition A^H* . 

From (10) and using (1), since E^ is (p^ x r^) and Y^ is (r^ x q) 
we have E = G = EA and Y = W = D®H* . The joining vectors are the 
columns of A which are the eigenvectors of . The columns of E 
(or equivalently G) are the spatial global EOF's and the rows of 
Y (or equivalently W) are the corresponding time coefficients or 
principal components . 


(12) 


I 

It is also shown in Appendix C that if (i=l , 2, . . . ,M) 

and k - k^ + . . . + kj^j are the smallest integers such that 


j=i 


and 




( 11 ) 


where «<' is a specified tolerance level, is the trace of , 

+ . . . + Tjyj , j^d^ are the eigenvalues of , and' 

the d. are the eigenvalues of YY' , then we can approximate V by 
J 

V at the 100 level, where 


V = g. w. + . . . + g_w_ 
^ ^ k k 


( 12 ) 


k 


( 13 ) 



3.0 ESMR DATA IN POLAR RESEARCH 


Passive microwave images of the polar regions obtained from 
the Nimbus-5 ESMR brightness temperature measurements have 
been shown to be a valuable source of polar information. The ESBflR 
receiver ( 1.55 cm) is generally useful in all-weather all-season 
situations, and, in particular, is able to detect the large con- 
trast between the brightness temperature of sea ice and the bright- 
ness temperature of open water. It is because of this contrast 
that the edge of the ice pack can be identified and information on 
sea ice concentration and ice type derived from these images. 

Scanning of the ESMR sensor provides complete spatial detail while 
continuous satellite coverage of these regions on a daily basis 
provides temporal detail (Zwally and Gloersen 1977). 

Gloersen et al (1978) characterized the time variation of the 
sea-ice concentration and multi-year ice fraction within the pack 
ice in the Arctic Basin through analysis of the ESMR microwave 
images and other data acquired using the NASA CV-990 airborne lab- 
oratory. The data were analyzed for four seasons during 1973~1975. 

These observations have shown significant variations in the sea-ice 
concentration in the spring, late fall and early winter. 

Carsey (I 98 O) used the ESMR microwave images in a study of the 

long-term and short-term behavior of the Weddell Polynya for the 

years 1973~1977. This polynya or ice-enclosed open area was observed 

during the 197^ » 1975* and 1976 winters. The behavior of this polynya 

margin and the regional ice concentration are interpreted in light 

of several oceanographic and meteorological theories explaining the 

circulation relevant to its origin, stability and role. He concludes 

that water column stability preconditioning alone is a necessary, but 

not sufficient condition for the existence of the polynya. 

___ ^ 

The Nimbus-5 spacecraft was launched 11 December 1972 into a nearly 

circular polar orbit (IO 89 x 1102 km) permitting complete global 

coverage every 12 hours. 

(4) The microwave radiation thermally emitted by an object is called its 
brightness temperature. It is expressed in units of temperature since 
the radiation emitted by a perfect emitter is proportional to its 
physical temperature for wavelengths in the microwave region. 

(Zwally and Gloersen 1977) 
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Crane et a.l (1982) used the ESMR data to determine the 
spatial and temporal patterns of change in microwave signa- 
tures of Arctic sea ice during a full annual cycle (1973/197^) • 
Interactions of ice conditions with the atmosphere are examined 
using grid point data for surface air temperature and atmos- 
pheric pressure. An EOF/principal component analysis is used to ex 
amine the major elements present in the microwave and atmospheric 
data. Principal components from these analyses are then used in 
a canonical correlation analysis to determine interassociations 
present hetween the ice and atmosphere in the Beaufort Sea and 
the European sectors on a synoptic time scale. 

Rayner and Howarth (1979) studied the areal extent and var- 
iability of Antarctic sea ice using the ESMR brightness temper- 
atures from 19 December 1972 to ^ June 1975 for a total of 219 
three-day time points with two data gaps (one from 26 February 
to 29 May 1973» "the other from 1 August 1973 to 5 September 1973)* 
The seasonal and interannual variations in Antarctic sea ice, 
e.g., the times of ice break-up and ice build-up, are important 
for fishing operations and coastal navigation. Also, these varia- 
tions influence global climate. 

In their method of analysis by Fourier decomposition Rayner 
and Howarth generate two series: one of maximum ice extent (min- 
imum latitude or MINI) and one of minimum ice extent (maximum 
latitude or MAXI) . Both these series with values at the integer 
longitudes follow the 155 K isotherm around the Antarctic contin- 
ent and are assumed to coincide with the 15?^ ice-sea isopleth 
which they take to define the ice-sea boundary (Zwally et al 1976) 
Land masses are initially masked out. The reason for using two 
series is that quite frequently, due to polynyas and large embay- 
ments within the ice pack, the 155 K isotherm crosses the same mer 
idian more than once. For a Fourier series representation of the 
ice-sea boundary the function must be single- valued as well as 
periodic . 

For the 219 time points both the MAXL and MINI series (lati- 
tude of the 155 K brightness temperature isotherm versus integer 
longitude) were expanded in a Fourier series each containing 180 



frequencies corresponding to the 3^0 longitudinal points. Results 
of their analysis reveal that the first seven harmonics are 
sufficient to account for more than '^0‘fo of the variance in the 
MINI boundary over the observation period and more than SOfo of 
the variance during the winter seasons. 

From their analysis of the amplitudes of the Fourier series 
with time they conclude that; (1) the time variation of the spa- 
tial mean latitude of the boundary is relatively smooth with a 
maximum occurring on 4 September 197^ and a minimum on 19 Feb- 
ruary 1975 (Figures 1 and 2) ; (2) the cold seasons are relatively 
long (205 days) and the summer season relatively short (l 60 days) ; 

(3) the time variation of the spatial mean latitude is asymmetric 
with the spring-summer removal of ice lasting 120 days and the 
autumn-winter increase lasting 150 days. 

Cavalieri and Parkinson (1980) have used Rayner and Howarth's 
KIINL series every ten degrees of longitude for the year 197^ and 
have compared the three-day averaged ESMR brightness temperatures 
with 1000 mb temperature and sea level pressure fields obtained 
from the Australian meteorological data set. In their work they 
also performed a Fourier decomposition - one on each of these 
three variables - and concluded that the first three harmonics 
were sufficient to account for most of the variance of the sea ice 
extent and temperature for any three-day period. Their results 
demonstrate an ice-atmosphere coupling of varying strength through- 
out the year. 

Zwally et al (I 98 I) have created a summary data set contain- 
ing Antarctic sea-ice conditions derived from the ESMR brightness 
temperature measurements for the years 1973 through 1976. The meas- 
urements have been mapped onto a polar stereographic grid enclosing 
the 50°S latitude circle. Sea ice concentrations have been calculated 
from the data for each grid element with an algorithm which uses an 
emissivity of 0.9 and an ice physical temperature estimate from 
climatological surface air temperatures. Monthly, multi-year monthly, 
and yearly maps of brightness temperatures and sea ice concentrations 
have been created. In their work they conclude that the microwave 
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brightness temperatures of Antarctic sea ice are predominantly 
characteristic of first-year ice with an emissivity of 0.92 at 
the 1.55 cm wavelength of the ESMR. 

A detailed analysis of Antarctic sea ice for the years 
1973 through 1976 as derived from satellite passive microwave 
observations has been prepared (Zwally et al I983) and is con- 
tained in an atlas which includes pseudo-color images of bright- 
ness temperature, sea ice concentration, changes in ice concen- 
tration, and multi-year average concentration. 
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3.1 EMISSIVITY. BRIGHTNESS TEMPERATURE. AND SEA ICE CONCENTRATION 


The Nimhus-5 ESMR receiver has a center frequency of 19 -35 
GHz (1.55 cm) with an IF bandpass of from 5 "to 125 MHz and is 
therefore sensitive to radiation from 19.225 to 19-^75 GHz, except 
for a 10 MHz gap in the center of the band (Wilheit 1972) . The 
radiometer is of the Dicke-type (Dicke 19^6) with a temperature 
sensitivity of 2°K and is fed by a phased-array antenna consisting 
of 103 waveguide elements each having an electrical ferrite phase 
shifter which step-scans across the sub-satellite track in 78 
positions. The total swath is + 50° from nadir (Wilheit 1972, 
Gloersen et al 1973) - Microwave images produced by the ESMR 
instrument have a spatial resolution of approximately 30 kilometers. 

Radiation (or brightness temperature) at the 1.55 cm wave- 
length is directly proportional to the received radiometric power 
since the Rayleigh-Jeans approximation holds. It is affected by 
high humidity and large water droplets in the atmosphere. Since the 
humidity is low and the number and size of water droplets are small 
in the south polar region, the received radiation can be assumed to 
be emitted from the earth's surface. Also, the radiation is not 
affected by darkness or cloud cover (Wilheit 1972, Wilheit et al 

1972). 

The basic equation for microwave radiometry is (Wilheit 1972, 
Zwally and Gloersen 1977) 


Tg = ^T ( 13 ) 

where Tg is brightness temperature, ^ is emissivity and T is the 
thermodynamic temperature of the emitting surface. 

For microwave radiation of wavelength 1.55 cm, the emissivity 
of sea water is approximately 0.4, first^year sea ice approximately 


(5) The Rayleigh-Jeans approximation for the intensity of thermal 
radiation from a blackbody works well at microwave frequencies 
(I -500 GHz) and at temperatures typical of the earth and its 
atmosphere (200 - 300°K) (Wilheit 1972). 
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0.95» and multi-year ice from 0.8 to 0.9 (Wilheit 1972; Gloersen 
et al 1973 » 197^) • More recent studies of first-year ice in the 
southern ocean have used a value of O .92 for first-year sea ice 
emissivity (Zwally et al I 98 I). 

For an area containing only first-year sea ice and open 
water, which is the case for most of the ice -influenced waters 
surrounding Antarctica, the brightness temperature Tg in a single 
instantaneous field of view (IFOV) is related to fractional ice 
concentration C (the fraction of the surface area within the IFOV 
which is covered by ice) by 


T 
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( 14 ) 


where € is first-year sea ice brightness temperature (assumed 
to be 235 K), ^ is open water brightness temperature (assumed 
to be 120 K), and A is the brightness temperature contributed by 
atmospheric water vapor (assumed to be 15 K over open water areas 
and 0 over ice) (Zwally et al 1976; Zwally and Gloersen 1977)* 
Neglecting variations of the ice physical temperature (T^)» sub- 
stituting the foregoing approximate values into (l4) and mult- 
iplying by 100 to express C in percent, we obtain 

C = (Tg - 135)?5 (15) 

Equation (15) is accurate to approximately 15?^ over most 
of the Antarctic region where sea ice predominantly has the 
microwave emissivity of first-year ice even though in some regions 
it may be older than first-year ice (Zwally et al I 983 ) • 
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4.0 ANALYSIS AND RESULTS 


In Section 4,1, characteristics of the data set are discussed, 
such as the period of observation, spatial coverage, temporal and 
spatial resolution, and missing data points. The minimum and maximum 
ice extents during the period of observation are shown, as well as 
time series plots of the spatial mean vector, its associated spatial 
standard deviation, and power spectra of these functions. 

In Section 4.2, the temporal mean and standard deviation ice 
brightness temperature maps for this time period are discussed. Char- 
acteristic features such as the Weddell Polynya, the edge of the pack 
ice, and the region occupied by ice year round, are identified. It is 
suggested that the advance and retreat of the ice edge over the annual 
cycle may be characterized by a conceptualy simple one-dimensional 
model with the gradient of mean ice concentration contours along a 
"growth/decay trajectory" on the mean map describing a nearly linear 
ramp function. 

An EOF decomposition of the sea ice brightness temperature anomaly 
field is discussed in Section 4.3. Gray scale/contour maps of the first 
10 E0F*s are shown, along with their corresponding time coefficients 
(or principal components) and variance spectra of the latter. In the 
interpretation of the B0F*s, the spatial pattern of anomalies over the 
entire grid of points can be considered to be a pattern of constructive 
and destructive interferences between the standing wave patterns in 
the EOF maps (Figures 6 through 10). In context with the one-dimensional 
model, the interference patterns between the EOF standing waves aerve, 
in part, to reconstruct the advance and retreat of the step-like ice 
edge through the annual cycle of growth and decay. A physical interpre- 
tation of the first four EOF's and their coefficients and variance 
spectra is given. 

In Section 4,4, the method of matrix partitioning in the time 
domain to approximate the "global" principal component structure is 
compared to the direct method of calculating global principal compon- 
ents. Calculations for the first three "pseudo" principal components 
are carried out for a tolerance level of 81^ (i.e., retaining 81^ of 
total sample variance) and are in good agreement with the first three 
global components, 
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4.1 DATA SET CHARACTERISTICS 


An extensive time series of Nimbus-5 ESMR images has been 
processed at Goddard Space Flight Center. ESMR brightness tem- 
peratures accumulated during each 3-day period since 30 Septem- 
ber 1973 were resampled to calculate 3~day average brightness 
temperatures at each of 293 x 293 grid points which are equally 
spaced when overlaid on the polar stereographic map projection 
illustrated in Figures 1 through 3» The average distance between 
grid points is approximately 30 km. The characteristics of this 
data set are more fully described by Zwally and Gloersen (1977 )• 

For the present analysis, an ensemble of ESMR sea ice Tg 
maps was selected for the time period extending from 30 September 
1973 through 25 May 1975* Of the 201 3-day time intervals spann- 
ing this 20 month sample period, 170 usable sea ice maps exist, 
and the 3I missing data points are distributed in relatively 
short time gapss the largest gap spans 6 time points, another 
spans 5 points, and the rest are smaller. The missing data in 
the time domain were treated using the Fourier analysis methods 
described by Murray (I98I). 

The space domain selected for analysis was restricted to 
the 15*597 oceanic grid points falling within the box labelled 
"Area Analyzed" in Figures 1 through 3* The raw data vector, at 
any time point t, thus contains the 3-day average brightness 
temperature at each of these 15*597 grid points (arranged in a 
fixed order). Recalling equation (15)* we may interpret each 
(Tg - 135) as ice concentration in percent at time t and the 
associated grid point. Missing data at individual grid points 
within any given map were filled by interpolating linearly over 
time at that location. 

Figures 1 and 2 illustrate the positions of the ice edge_, 
during this ensemble period, at the respective times of maximum 
and minimum sea ice extent (after Rayner and Howarth 1979)- Tg 
variations at open ocean grid points seaward of the maximum ice 

extent (Figure 1) are associated with atmospheric water vapor 
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and liquid water droplets in passing weather systems. These sig- 
nals typically range from 5 to 10 Kelvins. Figure 3 illustrates 
the geography of the. Antarctic Continent, major glaciers, con- 
tinental ice sheets, various islands, and South America, in re- 
•lation to the oceanic area covered by the ensemble grid. 

The fraction of the total analysis area covered by sea ice 
at any given time may be estimated approximately by averaging 
each data vector over the 15»397 grid points. The time series of 
the spatial mean vector and the associated time series of spatial 
standard deviations are plotted in Figure 4. Also illustrated in 
Figure k are the .least squares power spectra of these series 
(Vanicek 1969» 1971; We.lls and Vanicek 1978; Murray 1981). Both 
the spatial mean and spatial standard deviation are dominated by 
the annual cyc.le with a significant spectral peak associated with 
a half-year period. At the time of maximum ice extent, approx- 
imately 40^ to ^’5fo of the total oceanic area within the analysis 
area (Figure 3) is ice-covered. At the time of minimum ice extent, 
approximately to is ice-covered, using the approximation 
that (Tg - 135 ) is ice concentration in percent. 

4.2 TEMPORAL MEAN AND STANDARD DEVIATION OF SEA ICE BRIGHTNESS 
TEMPERATURE IN THE WEDDELL SEA (9/30/73 - 5/26/75) 

The temporal mean and standard deviation ice brightness 
temperature maps for this time series are Illustrated in the 
top two pane, Is of Figure 6. The 154 K mean brightness temper-^ 
ature isotherm and iso.line of 14 K standard deviation corres- 
pond c.lGBe.ly to the maximum ice extent (Figure 1). The open 
ocean region seaward of the maximum ice extent is character- 
ized by an average brightness temperature s.light.ly greater 
than 135 K, together with a standard deviation of approxi 
mate.ly 10 K due to water vapor and .liquid water in weather 
systems (Allison et al 1974; Wi.lheit et a.l 1975; Zwa.l.ly and 

Q.loersen 1977)- 
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The region occupied by sea ice year round, just to the 
east of the Antarctic Peninsula, is characterized by a high 
average brightness temperature and also a standard deviation 
of only 10°K. 

The 200 km x 1000 km area of low average brightness temp- 
erature near 0° longitude is a persistent, major polynya which 
has been studied by several investigators and called the "Weddell 
Polynya" (Zwally and Gloersen 1977; Rayner and Howarth 1979; Car- 
sey I 98 O; Parkinson 1983) • The ice pack completely covered the 
Weddell Polynya in the winter of 1973, but during the winters of 
197 ^ and 1975 this region was primarily open water with ice con- 
centrations rarely exceeding 1^% (bottom panel of Figure 11). 
Variation in ice concentration in this vicinity in 197^ was pri- 
marily due to a gradual change in shape of the polynya from Aug- 
ust through November 197^ (Carsey I 98 O). The standard deviations 
of grid points within the Weddell Polynya are corresponding low, 
being only slightly greater than those of the open sea beyond the 
maximum ice extent. 

The Weddell Polynya lies along the locus of the Antarctic 
Divergence, as it projects into the Weddell Sea. Enhanced diver- 
gent wind stress, due to inversion winds off the ice shelf, are 
thought to produce enhanced upwelling of warm deep water in this 
region and thus maintain the polynya as an open water area (Zwally 
et. al 1976 ; Zwally and Gloersen 1977). 

Parkinson (1983) » using mean climatological data as input, 
has successfully simulated a Weddell polynya through a sequence 
of numerical simulations of the ice cover, and concludes, "... 
the model obtains a full-scale, 0% concentration Weddell polynya 
which can be moved and prolonged by changing the wind forcing. 
Results seem to suggest that the Weddell polynya might arise in 
response to the winds, but then is maintained through oceanographic 
factors, including possibly a wind-driven divergence, or through 
the joint influence of the ocean and such feedbacks as the atmos- 
pheric warming which the polynya would be expected to induce." 
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In the mean field, the curved tongue-like protuberance 
emanating eastward from the Antarctic Peninsula represents the 
path of maximum growth and decay of the ice pack over the 
annual cycle. For purposes of discussion, the dashed line down 
the center of this tongue will be referred to as the ice 
"growth/decay trajectory". The gradient of mean ice concen- 
tration contours along the growth/decay trajectory may be char- 
acterized as a nearly linear ramp function. The ice growth/ 
decay trajectory curves around the Weddell Polynya, reflecting 
the low average ice concentration there, even though it closed 
over in 1973- 

The ramplike nature of the mean Tg contours along the 
growth/decay trajectory can be readily explained in terms of a 
simple one-dimensional model. Consider a horizontal axis from 
X = 0 to X = 1, along which ice concentration C(x,t) at any 
time t may be either 0 (open water) or 1 (ice). At time t = 0, 
C(0,0) = 1 and C(x,0) = 0 for x > 0. Over some period T, the 
ice edge is observed to advance from x = 0 (at time t = 0) at 
a uniform speed to arrive at x = 1 at time t = T/2, and then, 
to retreat at the same rate, returning to the initial con- 
ditions at time t = T. If we consider the mid-point at x = 0.5* 
C(0.5*t) is 0 for half of the cycle and 1 for the other half, 
with a mean (over T) of 0(0.5) = 0.5. Similar consideration of 
the mean at a few other points quickly reveals that C(x) is a 
linear ramp function sloping from C(0) = 1 to C(l) = 0. The 
mean concentration in this simple model is obviously similar 
to the contour gradient along the ice growth/decay trajectory 
in the mean ice concentration field (Figure 6 ) . 

If we now consider the standard deviation associated with 
the above one-dimensional model, Sq(x) = -x) , we obtain for 
a few selected points: S^(0) = S^(l) = 0, S^(.25) = S^(.?5) = 
0.^33* S^(.5) = 0.5. This is, therefore, a symmetric function 
with steep slopes near x = 0 and x = 1 , and a broad maximum at 
X = 0,5. Again, there are obvious similarities between the 
standard deviation curve of the model and the contours of stan- 
dard deviation in sea ice brightness temperatures along the 
growth/decay trajectory (upper right panel of Figure 6). 



From the above qualitative comparisons of geometric char- 
acteristics between the modeled and observed sea ice Tg means 
and standard deviations, we may conclude that much of the ob- 
served variance m.ay be explained by analogy to a simple 
"advancing step function" model of the ice pack. This inter- 
pretation is also supported by the time sequences of ice con- 
centration at selected individual pixels within the region in- 
fluenced by the ice pack (Figure 11 ) . In all three cases shown 
here, the transition from open water to ice-covered (and vice 
versa) occurs abruptly. Near the center of the pack, the en- 
semble average ice concentration - again approximated as 
(Tg - 135 )?^ - is close to 50%. And^at the pixel near the ice 
edge, the ensemble average ice concentration is only a few per- 
cent. When we consider the EOF decomposition of sea ice bright- 
ness temperature anomalies from the sample mean (Section ^.3 

below) , it will therefore be helpful to think in terms of re- 
constructing a progressive, slab-like advance and decay sequence 
through interference between "standing wave" EOF patterns, in 
combination with a linear mean slope extending from minimum to 
maximum ice extent. 


4.3 EOF DECOMPOSITION OF SEA ICE BRIGHTNESS TEMPERATURE ANOMALIES 
IN THE WEDDELL SEA (9/30/73 - 5/25/76) 

The usual objective of EOF analysis of geophysical obser- 
vations is to rotate the representation of the data to distances 
along a few p- dimensional directions containing most of the ob- 
served signal. This geometrical perspective is discussed in more 
detail here in Appendix A and elsewhere in many of the references 
cited above in Section 1 . If the original p-dimensional measured 
variables are truly random and uncorrelated with each other, then 
there will be no dominant dispersion (variance) directions and 
EOF analysis becomes a fruitless exercise. It is more frequently 
the case with geophysical observations, however, that strong corr- 
elations exist between the original variables (the measurement 
of these variables being made independently) , so that preferred 
variance directions do exist and can be illuminated through EOF 
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decompositions. Stated another way, EOF methods are often useful 
in situations where all p measurables respond to variations in 
some common forcing function, such as seasonal variation in solar 
radiative flux. The projections of the individual responses onto 
linear combinations of the original measurables (each such com- 
bination being a rotated direction) will influence the EOF coord- 
inates according to the importance of that process in the physical 
system being studied. 


In the case at hand, the original p dimensions are raw bright- 
ness temperatures at the 15 i 597 spatial grid points. A "direction" 
in this coordinate frame thus would be displayed graphically, 
either as a curve (if all grid points were indexed and plotted on a 
single abscissa with Tg as the ordinate) or as a two-dimensional 
gray scale/contour plot where the grid points are located on a map 
projection. The original measurement origin is thus located at zero 
brightness temperature at all grid points, a position far from the 
"centroid" of the sample, as given by the sample mean (upper left 
panel of Figure 6 ) . Before rotating to find principal axes of 
brightness temperature variations, it is therefore necessary to 
center the data matrix by subtracting the temporal mean from each 
observation. In geometric terms, this amounts to translating the 
origin to the position of the sample mean or centroid of the data. 
The resulting data matrix V in (1) thus contains column vectors of 
brightness temperature anomalies, v., relative to the temporal mean^ 
for 170 time points. These anomaly values, in Kelvins, are now di- 
rectly interpretable as ice concentration anomalies in percent. 


The EOF's defining the principal directions of anomaly var- 
iations in ice concentration are the eigenvectors of the spatial 
sample covariance matrix 


S = vv'/(i-i) 


VV/169 


where q = I 70 is the sample size in the time domain. 

Since the mean is subtracted from each observation there are 
at most 169 degrees of freedom and, therefore, the rank of S, or 
the number of nonzero eigenvalues, is less than or equal to I69. 
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The elements along the main diagonal of S are the sample 
variance of brightness temperature measurements at each assoc- 
iated spatial grid point. Therefore, the square roots of these 
elements yield the sample standard deviation shown as a gray 
scale/contour map in the upper right panel of Figure 6 and dis- 
cussed in Section 4.2. The trace of (sum of all the elements 
along the main diagonal) is total sample variance, a quan- 
tity which is invariant under the rigid axis rotations assoc- 
iated with a complete EOF decomposition (Section 2 and Appen- 
dix A) . 

In order to reduce computation time and storage require- 

A 

ments, the eigenvalues and eigenvectors of S were calculated 

by first determining the eigenstructure of the smaller time 
domain scatter matrix V'V which has the same eigenvalues as W . 
The desired spatial eigenvectors and temporal principal compon- 
ents were then obtained as explained in Section 2 and Appendix B. 

A determination that only the first k eigenvalues are 
statistically greater than zero amounts to a statement that the 
first k principal components represent signal^ and the remainder 
represent noise. Preisendorf er et al (1982) have recently reviewed 
this topic and provide a systematic set of rules for selecting 
significant eigenvalues, but ambiguities remain even so. 

One family of eigenvalue selection rules is based on slope 
variations in a plot of the logarithm of eigenvalue (LEV) ver- 
sus eigenvalue number. Craddock and Flood (19^9) and Craddock and 
Flintoff ( 1970 ) » fon example, truncate the "signal" eigenvalues at 
the point beyond which the LEV curve may be approximated by a 
straight line. Other criteria for determining truncation level in- 
clude retention of eigenvalues exceeding of the total sample 
variance (Kaiser i 960 ), retaining only eigenvalues which are "much 
larger than others" (Beale et al 196?) , retention of eigenvalues 
which cumulatively account for 99 ?^ of the total sample variance 
(Mueller 1976), and subjective determinations based on inspection 
of the LEV diagram (Rinne and Jarvenoja 1979) • In all cases, some 
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element of subjective interpretation is involved in the trunc- 
ation decision. 

A 

The LEV diagram of the sea ice covariance matrix S is 
plotted in Figure 5. By the Craddock et al criterion, the first 
85 EOF's should be retained to provide a "complete" decomp- 
osition of the sea ice anomaly field. This truncation level 
would retain .^2fo of total sample variance. The first 10 
eigenvalues of S and their respective and cumulative contribu- 
tions to the total sample variance are compiled in Table 1. Each 
eigenvalue equates to variance in the principal component (dis- 
tance) in the direction defined by the associated eigenvector 
(EOF). Nearly ^2% of total sample variance in ice concentration 
is resolved by a principal component representation in these 10 
directions. By using either Craddock's criterion or Mueller's 
we would retain the first 85 EOF's. By neglecting those eigen- 
values which individually account for less than 1% of total var- 
iance we would retain only the first 7 EOF's. 

The issue of truncation would become critical in a study 
which used a trxincated EOF decomposition of the data in con- 
junction with a numerical sea ice model, or in the calculation 
of cross-covariance spectra between truncated representations of 
sea ice concentration and another variable, such as atmospheric 
pressure^ Truncation is less critical in purely descrip- 

tive analyses of organized structure in space/time variability 
associated with a geophysical field. The present analysis is arb- 
itrarily confined to a descriptive presentation and interpretation 
of the first 10 EOF's and principal components of sea ice concen- 
tration anomalies in the Weddell Sea. Furthermore, detailed dis- 
cussions focus primarily on only the first ^ EOF' s^ which are 
assumed to represent the dominant normal modes of variation in the 
field. 

The first 10 EOF's are illustrated as gray- s cal e/c on- 
tour maps. in Figures 6 through 10. Shown beside each EOF map are 
the corresponding principal component and the variance spectrum of 
the principal component time series. In all EOF maps, major anom- 
aly amplitudes are confined to grid points falling between the 

maximum ice extent boundary and the shoreline of Antarctica. In 

-28- 



many of the EOF's, very small amplitude anomaly patterns occur 
over the open ocean areas beyond the maximum ice boundary. These 
anomalies are due to water vapor and liquid water fluctuations 
along the storm tracks through this region during this ensemble 
period . 

The k th principal component at any time t represents the 
amplitude of the normalized (unit length) spatial anomaly pattern 
illustrated in the k th EOF map. The interplay between an EOF and 
its time varying principal component can be most easily under- 
stood through a specific example. Examine EOF/principal component 
#2 (Figure 6) and note that the EOF pattern indicates a negative 
anomaly of < -0.016 at position (60°S, 15°E) and a positive anom- 
aly of +0.012 at position (65°S, 40°W) . On day 250f the 2 nd prin- 
cipal component has an approximate value of +2000 K, which yields 
2 nd EOF anomaly contributions of <-32 K at (60°S, 15°E) and 
+ 24 K at (65°S, 40°W) . Conversely , on day 400 the 2 nd principal 
component has an approximate value of -1000 K, yielding 2 nd EOF 
contributions of > l6 K at (60°S, 15°E) and -12 K at (65°S, 40°W) . 
Reconstruction of the time modulation of the overall brightness 
temperature anomaly at these two locations requires, of course, 
that the contribution of all significant EOF’s be summed. Thus, 
the spatial pattern of anomalies over the entire grid can be con- 
sidered to be a pattern of constructive and destructive interfer- 
ences between the standing wave patterns in the EOF maps(Figs 6-10). 

Recall now the discussion of the simple one-dimensional 
"advancing step function" model introduced in Section 4.2 above. 

In the context discussed there, the interference patterns between 
EOF standing waves serve, in -part, to reconstruct the advance and 
retreat of the step-like ice edge through the domain influenced 
by growth and decay of sea ice through the annual cycle. The 
abrupt transition from water to ice (and the reverse) observed at 
individual pixels (Figure 11) is consistent with this interpretation 

In general, the 1^ EOF and principal component describe a 
long wavelength modulation of positive and negative anomalies, foll- 
owing a spatial pattern which is qualitatively similar to the shape 
of the standard deviation of the one-dimensional model (Section 4.2) 

Near the time of maximum ice extent (Rayner and Howarth 1979), the 
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1st EOF anomaly monotonically increases Tg relative to the mean, 
with maximum contributions near the center of the ice pack. Con- 
versely, near the time of minimum ice extent, the 1st EOF anom- 
alies are monotonically negative and tend to reduce Tg over the 
ice-influenced region towards the uniformly low values associated 
with open water. 

The 1st principal component varies smoothly between times of 
minimum and maximum sea ice coverage in concert with previously 
published descriptions of the annual cycle of ice growth and decay 
in the Weddell Sea (Rayner and Howarth 1979: Zwally and Gloersen 
1977). Not surprisingly, therefore, the variance spectrum of the 
1 st principal component is dominated by the annual cycle, with 
smaller spectral peaks associated with a half-year period and a 
120 day period. (The occurrence of peak power at 3^8 days, rather 
than 385 days, is almost certainly an artifact of the ensemble 
period being only 20 months. This series is too short to allow 
accurate estimates of a variance spectrum dominated by a 12-month 
cycle . ) 

The 2 nd principal component time series is phase shifted 
from the lst_ in such a way that interferences between l£t and 2 nd 
EOF anomaly patterns will improve the representation of the advance 
and retreat of the ice edge during periods of growth and decay 
(Figure 6) . A similar conclusion may be drawn from inspection of 
interactions between the first EOF/principal component and numbers 
3 and 5 through 10. In each case the largest principal component 
amplitudes (+ or -) occur during periods of ice growth and/or 
decay. Also, in each case the dominant spatial pattern in the 
higher EOF may be seen to produce interferences with the pattern of 
the mean and 1st EOF which reproduce a step-like manifestation of 
the ice edge. 

Many of the strongest anomaly patterns occur in regions where 
the ice pack undergoes extremely rapid growth and decay during the 
transition phases. This can be seen by comparing the EOF anomalies, 
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in this light, with 197^ monthly average sea ice extent contours 
published by Cavalieri and "ferkinson (1981) and included here as 
Figure 12. During the ice pack growth period in March through May, 
the spacing between monthly contour intervals is large in the 
western and central areas of the Weddell Sea (Figure 12a), Compare 
this especially with the anomaly patterns in the 2 nd and 3 rd EOF's 
and the rapid positive growth in their principal components during 
this period (Figures 6 and 7). Then^ during the early retreat of 
the ice pack in November and December, a large contour interval 
outlines the rapid expansion of the polynya, followed by an abrupt 
retreat to nearly minimum extent in January 1975 (Figure 12b) . The 
breakup of the pack ice followed a more complex sequence of spatial 
patterns than did the ice growth, and is manifested by rapid posi- 
tive to negative excursions in most of the principal components 
(Figures 6 - 10) . 

EOF #4 and its principal component time series (Figure 7) 
uniquely represent interannual variability in the Weddell Polynya 
between 1973 and 1974. In the winter of 1973 » when the polynya was 
absent until early in the spring breakup of the ice pack, positive 
brightness temperature anomalies prevailed here. Then, in 1974, the 
polynya remained open throughout the winter, creating a negative 
Tg anomaly relative to the 20-month mean. The anomaly amplitudes in 
EOF #4 also indicate a more subtle effect associated with the band 
of -0.012 anomaly extending zonally eastward from the tip of the 
Antarctic Peninsula (Figures 3 and 7)* In 1973* when the 4 th prin- 
cipal component had a large positive value (polynya absent), the 
maximum ice extent boundary in these longitudes was displaced anom- 
alously to the south. In 1974, when the 4 th principal component had 
a large negative value (polynya present), the maximum ice extent 
boundary was displaced anomalously to the north. A much longer time 
series must be analysed to determine whether there is significant 
evidence of an association between presence of the Weddell Polynya 
and maximum ice extent. In the present sample this apparent assoc- 
iation could be a purely coincidental contrast between two winters. 

The primary roles of EOF's 2,3, and 5 through 10 in reconstruct 
ing the movement of the ice edge during growth and decay phases is 
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also manifested in the time histories of each principal com- 
ponent's percentage contribution to the spatially averaged mean 
squared anomaly (equation (8) and Figures 13 and l4). The 1 st 
EOF dominates the mean squared anomaly during periods of both 
maximum and minimum ice extent and falls abruptly during growth 
and decay phases. At the times of ice minima, contributions 
from BOF's 2 through 10 are negligible. During the periods of 
maximum ice cover EOF #4 is a significant contributor, repre- 
senting the markedly different manifestations of the Weddell 
Polynya in 1973 and 197^. The situation is quite different dur- 
ing the ice growth and decay phases. The mean squared contri- 
bution of principal component #1 falls steeply to zero, and then 
rises just as steeply, as the sign of the principal component 
reverses. Coincident with these "notch-like" signatures, the 
mean squared anomaly signals associated with principal compon- 
ents 2,3» and 5-10 appear as sharp pulses of energy. During 
each such pulse, the shorter wavelength anomaly patterns of the 
associated EOF's make strong contributions to the reproduction 
of the ice edge. 

The variance spectra of the first 10 principal components 

are all dominated by peaks associated with periods which are a 

significant fraction of the approximate 600 day ensemble period. 

As previously mentioned, this time series is much too short to 

attribute a high degree of accuracy to the spectra illustrated 

in Figures 6 through 10. Nevertheless, the dominant peaks can be 

clearly associated with well-known processes, including the 

“1 

365 day ( 0.003 day ) annual cycle , (particularly in spectra of 

_ A 

principal components 1 and 2), a I 60 day (O.OO 6 day ) warm 
season and 205 day ( 0.005 day ) cold season (Zwally and Gloersen 
1977; Rayner and Howarth 1979) 1 and a 120 day (O.OO 8 day”^) per- 
iod of ice removal (Rayner and Howarth 1979) • 
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4.4 EOF COMPUTATIONS USING PARTITIONS IN THE TIME DOMAIN 


It was noted in Section 4.3 that the actual EOF computa- 
tions were done using the smaller time-domain scatter matrix 
V'V» and the results transformed for application to the spatial 
covariance matrix ^ using the relationships explained in Sec- 
tion 2 and Appendix B. The 170 x 1?0 matrix V V is small enough 
to allow practical direct computation of its eigenvalues and 
eigenvectors using standard computer subroutines. Were sample 
size to grow much larger, however, direct eigenvalue computations 
would become impractical and the techniques outlined for separate 
computations of eigenvectors and eigenvalues of smaller partit- 
ions, and subsequent joining of the results (Appendices C and D) 
would become essential. Therefore, this data set was used to test 
the accuracy of global principal components estimated by joining 
truncated sets of principal components associated with the smaller 
dimensioned scatter matrices. 

The data matrix (V) was partitioned in the time domain into 
5 sub-matrices, each of dimension 15*597 x 3^ and the eigenstruc- 
ture of each of the 3^ x 3^ ~ l*-*»5) scatter matrices 

determined. Using the joining procedure outlined in Appendices C 
and D and retaining 8lfo of the total sample variance, a truncated 
set of 9 principal components (and associated EOF's) from all 
these partitions was joined to approximate the global principal 
component structure. These approximations (or pseudo -principal com 
ponents) were then compared with the global principal components 
calculated directly without truncation. The first three pseudo- 
principal components are plotted as open circles in Figure 6 and 7 
The 1^ pseudo -principal component series is virtually indisting- 
uishable from the directly computed one, and the comparisons are 
nearly as good for the 2 nd and 3 rd principal components. Clearly, 
this approach will be usable, given proper retention of total 
sample variance in the truncated computations, for obtaining acc- 
urate estimates of at least the lower ordered EOF's from indefin- 
itely large ensembles of sea ice brightness temperatures. 


- 33 - 



A similar computation was carried out by partitioning 
the spatial grid into three longitudinal sectors. Equally 
good approximations to the global EOF's were obtained with 
this approach to partitioning of the sample. The results from 
the spatial partition computations allow some interesting 
comparisons to be made concerning the global importance of 
locally important anomaly patterns, particularly in the higher 
ordered EOF's. However, that type of interpretation is beyond 
the scope of the present paper, and the results of spatial 
partition computations will not be presented here. 
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5.0 DISCUSSION AND CONCLUSIONS 


A matrix partitioning scheme has been presented for approximat- 
ing the eigenstructure of a large sample covariance matrix. The data 
array, a field of measured anomalies (of some physical variable) 
relative to their time averages, may be partitioned in either the 
time domain or the space domain. Eigenvectors of the smaller dimen- 
sioned covariance matrices associated with the partitioned data 
sets (and their principal components) may be calculated independent- 
ly and joined to approximate the eigenstructure of the larger co- 
variance matrix associated with the unpartitioned data set. The 
accuracy of the desired approximation (in terms of retaining a giv- 
en percent of total sample variance) and the magnitudes of the 
largest eigenvalues from the partition covariance matrices deter- 
mine the number of local EOF's (eigenvectors) and principal com- 
ponents to be used in the joining process. 

This method is shown to allow accurate estimation of spatial 
EOF's for time series of satellite image data where there is a 
large number of spatial grid points. The spatial covariance matrix 
associated with any ensemble of satellite images is far too large 
to allow direct computation of its eigenvalues and eigenvectors, 
the latter being the desired spatial EOF's. In many cases, however, 
the number of images in an ensemble (its time dimension) is suffic- 
iently small to permit direct computation of the eigenvalues and 
eigenvectors of the temporal scatter matrix. The latter eigenvec- 
tors can be easily transformed to obtain the eigenvectors of the 
spatial covariance matrix. As a result, great computational sav- 
ings are realized. 

There are a number of satellite data ensembles where the 
dimensions of both the time and space scatter matrices are too 
large to permit direct computation of the eigenstructure of either 
matrix. In such cases, the above mentioned matrix partitioning 
scheme must be employed to break the computational problem into 
manageable pieces. The methods by which the partition results are 
joined are described in detail in Appendices C and D. 

The accuracy and efficiency of the partitioning approach were 
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tested using a 20-month ensemble of Nimbus-5 ESMR measurements of 
sea ice brightness temperature emissions in the Weddell Sea. The 
time series of sea ice maps extended from 30 September 1973 
through 25 May 1975 and included I 70 3~day composite average maps 
(leaving 3 I reasonably well distributed time gaps). The spatial 
representation was restricted to 15i597 oceanic grid points in the 
Weddell Sea and adjoining Antarctic Circumpolar Current regime. 
Spatial EOF's were first computed directly, using the I 70 x I 70 
time scatter matrix as a computational expedient. Then, the first 
3 EOF's were computed indirectly using partitioning in the time 
domain, and in another case, the space domain. Both partitioning 
methods yielded comparable results and the calculated EOF's and 
principal components agree remarkably well with those computed 
directly. It is concluded that this method is extremely well suited 
for analysis of data sets much larger than the one examined here. 

The first 10 spatial EOF's, together, explain 92fo of total 
sample variance in sea ice brightness temperature (or percent ice 
concentration) anomalies relative to the ensemble mean. The mean 
and standard deviation are shown, by qualitative comparison with a 
simple 1-dimensional model, to be consistent with those of a con- 
ceptual model of the ice pack as a 2-dimensional, quasi-uniform 
slab which grows and shrinks in spatial extent over the annual 
cycle of ice growth and decay. In this context, the dominant anom- 
aly patterns in the spatial EOF's and the phase relationships of 
the principal components may be interpreted as a standing wave de- 
composition of the advance and retreat of the ice edge. At any 
time between maximum and minimum extent, the ice edge is represen- 
ted by a linear pattern of interference between the ensemble mean 
ice concentration and the standing wave anomaly patterns in the 
EOF fields as modulated in amplitude by the time-varying principal 
components . 

The fourth EOF/principal component uniquely accounts for the 
interannual difference in the Weddell Polynya, which was absent 
during the winter of 1973 arid present throughout 197^- This inter- 
annual contrast accounts for of total sample variance which is 
a highly significant fraction. Recall that the full variance assoc- 
iated with the annual cycle is present in the covariance matrix 
analyzed here! The spatial anomaly pattern in the fourth EOF also 
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shows that the presence of the polynya in 197^ was accompanied 
by an anomalous northward protrusion of the maximum ice extent, 
and the reverse in 1973 when the polynya was absent. A much 
longer time series would be needed to determine whether this 
apparent association of the polynya with maximum ice extent 
east of the Antarctic Peninsula occurs regularly, or was merely 
a coincidental artifact of the two years examined. 

Power spectral analysis of the principal components reveal 
periods which can be related to the seasonal cycle of sea ice 
growth and decay in the Weddell Sea, harmonics of this cycle, the 
cold season (205 days), the warm season (l60 days), and the dur- 
ation of spring-summer ice removal (120 days) as reported in the 
literature. 

The first four EOF' s and their components can be considered 
the dominant or principal modes of variation in the ice field, 
accounting for of the total information content in the data 
(field variance) . 
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APPENDIX A 


GEOMETRICAL INTERPRETATION OF PRINCIPAL COMPONENTS 


Principal components can be viewed geometrically as the pro- 
jection of q points in a p-dimensional Euclidean space onto k 
orthogonal vectors of a "best fitting" k-dimensional subspace 
(k-p) . Equivalently, they are linear functions of q measurements 
on p variables which possess certain optimal properties. The sub- 
space is "best fitting" (Pearson 1901) in the sense that the sum 
of squares of the perpendicular distances from the points to the 
subspace is a minimum. The motivation for projecting the points 
onto a subspace of lower dimension is to represent the q measure- 
ments by less than p linear functions in such a way that there is 
a minimal loss of information content. The minimization procedure 
is equivalent to determining k orthogonal axes or directions such 
that there is maximum dispersion along the first principal axis, 
a second largest dispersion along the second principal axis, etc. 

The dispersion is a minimum along the k th principal axis (the dis- 
persion is a measure of the variability and is the sum of squares 
of the perpendicular distances of the points to a line passing 
through the center of gravity of the points) . The solution to the 
problem of finding a "best fitting" subspace to a scatter of points 
in a higher dimensional space involves the eigenvalues and eigen- 
vectors of a matrix and was first investigated by Pearson (I 90 I) 
and Frisch (1929) . 

As a simple example, consider the scatter of points in Figure 
A-1. The " best fitting" one-dimensional subspace passing through 
the center of gravity C of the 7 two-dimensional points is the line 
LpQ (obtained by minimizing the sum of squares of the perpendiculars 
from the points to the line) . The " standard" least squares straight 



^ ' The development here follows that of Rao (1964) who gives an ex- 
cellent treatment on the use and interpretation of principal com- 
ponents in applied research and who discusses many of their op- 
timal properties. Kshsirsagar (1972) is also an excellent refer- 
ence on the subject and the proof of one of the optimal properties 
found in this appendix is due to him. The appendix is self-contained. 



line fit to these points is the line (obtained by minimizing 

the sum of squares of the vertical distances from the points to 
the line). Note that the dispersion of the points is a maximum 
in the direction of Lp^ (i.e., the second moment about a line 
passing through the center of gravity is the greatest in this 
particular direction). The second principal axis (not shovm) is 
orthogonal to Lp^ and is in the direction of Ljg. , also passing 
through C . 


y 


Center of 
Gravity 


T, 


First principal axis 

Perpendicular from 
point to subspace 


Least squares 
straight line L-j.^ 


. A i V 

1 Vertical from 

point to line 




Figure A-1 


We can represent q points in a p-dimensional Euclidean space 
by the matrix X 


X = 


^11 • * * ^Iq 


pi 


• • • 


X 


PI 


(Xp I ... I Xq^) ^^i j ^ 


(A.l) 


where the point, x., is the column of X. 

J 

The center of gravity of the q points is x where the i th 
coordinate is 


Xi = (1/q) 


3=1 


(i l,2,.«.,p) 


(A. 2) 





The points measured from their center of gravity, 

= (x.. - X. ) (i = 1,2 j = 1,2, ...,q) can be written 

1 J 1 


V = 


^11 • * * ^Iq 


^pl • • • V 


= <^1 


(A. 3) 


The total dispersion (or scatter) matrix S is 


S = VV = (s. .) 

X J 


(A. 4) 


and is positive semi-definite with p eigenvalues d^-dg- • • .-dp-0 
which are the roots of 


S - dl =0 
P 


(A. 5) 


Corresponding to each d. there is an orthonormal eigenvector e. 

J J 

which satisfies 




(A. 6) 


or in matrix notation 


E'SE = D = diag(d^...d ) 


(A. 7) 


where 


= TT« 


(D 

(D 

(A. 8) 

E'E = I 

P 

(A. 9) 



and 


I is the identity matrix of order p. 
P 


Equation (A.?) can be written 


S = d^e^^e^ + 


d^e^el 
P P P 


(A. 10) 


■•■••• + ®p®p (A. 11) 

The representation of S and in (A. 10) and (A.ll) is known 
as the spectral decomposition of the matrix S and respectively 
(Good 1969 and Appendix B) , 

By definition a k-dimensional subspace is specified by a point 
(or origin) and k orthogonal vectors. As previously ihentioned Pear- 
son defined a "best fitting" subspace as one in which the sum of 
squares of the perpendicular distances from the points to the sub- 
space is a minimum. Such a subspace passes through the center of 
gravity of the points. 

Let h^ , h^i . . . I hj^ be k orthogonal vectors which together with x 
specify the k-dimensional subspace. Without loss of generality we may 
assume the h^ are orthonormal vectors. The square of the perpendicular 
distance from the i^ point to the subspace is then 

1? = (x^-x) • (x^-x) - (|h^ (Xj^-x)] ^ + . . . 

(A. 12) 

Summing (A. 12) over all q points we obtain 



trace(S) - trace(H'SH) 


(A. 13 ) 


where 


H = (h^...h^) = (h..) 


(A. 14) 


H-H = 


(A. 15) 


Since the eigenvectors e . are orthogonal they are linearly 

3 

independent and form a basis for the p-dimensional space. Therefore 
we can write h. as a linear function of the e. 


^i " c^^e^ + . , . 


-f + c- e 

ip p 


(i ~ l»2j...|k) 


(A. 16) 


or in matrix notation 


H = EC 


where E is given by (A. 8) and 


C = 


All * 

\°kl • 



(A.l?) 


(A. 18) 


with 

CO’ = 


(A. 19) 


Using (A. 7) I (A. 9). and (A.l?) trace(H'SH) becomes 


trace (H'SH) 


trace(CDC ' ) 



• • • 


• • • 


(A. 20) 


di(cii + 


+ c^i) + 


+ d ( + 
p' ip 




where each of the coefficients of is - 1 and the sum of all the 
•coefficients is k by (A. 19)* 

In order to maximize trace (H'SH) the optimum choice for the 
coefficients is unity for d^, ...,dj^ and zero for the others. This 
means c.. = 1 for i = l,2,...,k, c. . = 0 for all other i and j 

^ ^ tJ 

which is possible by choosing h^ = e^ for i '= 1,2 k. 

Therefore 


max trace(CDC') 


d^ "i* i i • 


(A. 21) 


and the minimum of Q in (A. 13) is 


min Q = 
^i 


^k+1 ^ ^p 


(A. 22) 


which is attained when h^ = e^ (i = l,2,...,k) 


The actual representation of points in the "best fitting" 
lower dimensional subspace is, letting = (e^,...,ej^) 


Kh ■ ■ ■ 


= 


• • • 


(A. 23) 


e,'x. . . . elx 




k^q 


The first row of E^X gives the best one-dimensional representation 
of the points} the first two rows give the best two-dimensional rep- 
resentation, etc. 


In similar manner 


Y = Ej^V = 






(A. 2^) 


gives the "best fitting lower dimensional representation of the 
deviations of the points from their center of gravity. Here 
(i = l, 2 ,...,k) is the i^ principal component as used in the 
text . 


Another optimal property of the first k eigenvectors of S 
should he noted. 

In the sense of least squares the sum of the first k terms 
in the spectral decomposition (or singular value decomposition) 
of S gives the hest matrix of rank k that approximates S . 

Thus, if B is any matrix of rank k then 


min jjs ~ B jj = 
B 

^^k+l ^p^^ 

(A. 25) 

where the Euclidean norm of A, j|Aj.| 

, is 

D 


IN = 

‘ 

(A. 26 ) 


See Appendix B. 



APPENDIX B 


RELATIONSHIP BETWEEN SINGULAR DECOMPOSITION 
AND EOF/principal component analysis 


Let V Le a real (p x q) data matrix of rank r. Using the 
singular decomposition theorem (Good 19^9) we can write 


V = 


+ e d^f 
r r r 


(B.l) 


where the column vectors e- and f • are respectively the left 

J J 

singular vectors of V ( or eigenvectors of W ) and the ri^ht 
singular vectors of V ( or eigenvectors of V'V ), and the d^ 
are the singular values of V or the square roots of the positive 
eigenvalues of VV ( or of V'V ) with 


d! - dl ^ 


^ d^ 


> 0 


(B.2) 


and 


e.' e . = f f . = X- . 

J 1 J ij 


(B.3) 


If p is the number of spatial points and q the number 
of time points and the row means ( time averages ) have been 
removed from the data array V, then the e . are sometimes re- 
f erred to as spatial EOF's ( empirical orthogonal functions ) 
and the f. as time EOF's ( Lorenz 195^ ). 

J 

Using (B.l) and the orthonormal property of e. and f . in 

J J 

(B.3) we can write 

V = f^dfe' + ... + (B.4) 

VV = e^d^e^ + ... + e^d^e^ (B.5) 

V'V= + ... + qd^V (B.6) 

To show that e. is an eigenvector of VV corresponding to d^ 

J J 



we use (B. 3 ) ^nd (B.5) 


(vV)e. = d.ej 


To show that fj is an eigenvector of v' V corresponding to 
d. we use (B.3) and (B.6) 

J 




Define the principal component y- corresponding to e. as 

u J 


y ■ = e . V 


Then using (B.l) and (B. 3 ) becomes 

J 


y . = d®f 
J J J 


and we see that y'. is also an eigenvector of v'v corresponding 

J 

to d. (not normalized). 

J 


In similar manner define the principal component z • corres 

J 

ponding to f • as 
J 


^3 = 


Then using (B.3) and (B.4) z. becomes 

J 


“3 “ '^3®3 


(B.7) 


(B.8) 


(B.9) 


(B.IO) 


(B.ll) 


(B.12) 



and we see that z'- is also an eigenvector of VV’ corresponding 

3 

to d. (not normalized). 

Using (B.3)i (B.IO), and (B.12) we can write 


V = z’. = d, J, , 


V • V . = z . z 


3 1 13 


and we see that the y. and z. row vectors are also orthogonal 

3 3 


Using (B.IO) and (B.12) in (B.l), (B.4), (B.5), and (B.6) 
we can write the singular decomposition of these matrices in terms 
of the eigenvectors of VV' and V' V and their principal components. 



V = 

e^y^ + . . . + e^^y^ 

(B.14) 


V' = 

fl^l ^r^r 

(B.15) 


vv' = 

z^ z^ + . . . + z^z^ 

(B.16) 


v'v = 

yjy ^ + . . . + y^y^ 

(B.l?) 

The trace of 

vv' (or 

of V'V) is 


Trace (VV* ) 

= Trace(yjy^) + ... + Trace(y^y^) 



= Ill'll 

2 + ... + |yrf 



= yisi 

+ . . . + yr^r 



- d^ + d„ + . . . + d 
12 r 


(B.18) 



where jjx jj is the Euclidean norm of x. 

We can use the first k - r terms in the singular decomp- 
osition of the matrices in (B.l), (B.4) , (B.5)> and (B.6) to 
approximate them. Thus, for V we have the approximation V 


V = V 


i 1 . 

p ^ * + 4- p * 


= e^y.j_ + . . . + e^^y 




(B.19) 


In the sense of least squares V gives the best matrix of 
rank k that approximates V. That is, (V - V) has minimum norm. 

In order to see this, since the pq matrices e.f. (i = 1,2 p: 

j = 1,2 q) form a basis for the vector space of (p x q) 

matrices, we can write an arbitrary approximation for V with 
rank k, or as 


a 




1. j 


o^. .e.f. 

ij 1 J 


(B.20) 


as 


Writing the complete singular d-ecomposition of V (Good 19^9) 


V 



i. j 


y -e- f *4 

'ij 1 J 


we have 



r 

E (4Xi.i - 


i. j=l 


(B.21) 


(B.22) 


and the optimum choice for 



is 



lU. . 

1 .1 


(B.23) 


(i 



This gives V„ = V as the best matrix of rank k which 

Ci 

approximates V. The minimum attained is 


||v - vll^ = + ... + d 


k+1 


(B.24) 


Since the sum of all the d. (j = 1,2 r) is the trace T 

of W , each term in the expansion of V, e.g., the explains 


dj/(di + . . . + 100 % 


of T. 


In similar fashion it can be shown that the first k - r 
terms in the expansion of S = W in (B.5) gives us the best 
matrix in the least squares sense which approximates S. That is, 
if B is any matrix of rank k then (S - B) has minimum norm and 
equation (A. 25) in Appendix A holds. 

The above illustrates the dual role of eigenvectors and 
principal components in the V and V spaces. If e- is an ortho- 
normal eigenvector of VV corresponding to the eigenvalue d^ 0 
with principal component y. = e'-Y, then f . = d^^y'- is an ortho- 
normal eigenvector of VV. If f • is an orthonormal eigenvector of 

J 

VV corresponding to d - ^ 0 with principal component z- = f '• V , 

J ‘ J J 

then e . = is an orthonormal eigenvector of VV . 

J J J 



APPENDIX C 


RETATTONSHIP BETWEEN A PARTITIONED AND A NON- PARTITIONED MATRIX V 

Let V be an arbitrary real data matrix with dimension (p x q) 
and rank r - q - p. Then we can write (Appendix B) 


V = ED^F' = EY 


+ ... + e^d=f' 




+ e y 
rJ'r 


(C.l) 


with 


E'E = F’F = I^ 

E = (e^...e^) 

F = 

D^ = diag(dj. . . dj) 

YY' = D = diag(d^ . . .d^) 

where the e- and f. (j = l,2,...,r) are respectively the left and 
tJ iJ Jx. 

right singu.lar vectors of V and the d^ are the singular values of V, 

J 

Since the rank of V is r 




• • • 


0 



Partition V into M matrices with dimension x q). and 
rank r^ ^ q where + . . . + = p, r^ + . . . + r - r, and use 
the singular value decomposition for each ( i = ) 


V 


lv,\ 


^Vm/ 


El 


= e. y. + . . . + e_y_ 

i 1 . P P 


(C. 


where r = r^ + . . . + rjyj and 


^i 


E. D?F! 
11 1 


E. Y. 

1 1 


= ie^(iyi) + ... + le^ (iy ) 


(i = 1,2, , M) 



The e^ vectors (columns of E) form an orthogonal set. Let 
Pj_ = p^_^ + p^, Pq = 0, and r^ = r^_^ + r^, = 0. Then we can 
write 




i®kj 


for 1 - j - and 

Pi-i< k - Pi 
1 - 1 M 
0 otherwise 


where j^e^. 
le. of E. 


is the k th element of . e • 
are the eigenvectors 


and the column vectors 
from the i;^ subdivision. 


Using the orthonormal property of the 


. e . 

1 J 


vectors 





The y • vectors (rows of Y) are related to the -y • (principal 

J cl 

components corresponding to the ^e^ or row vectors of Y^) 




l^j (1 j - Pi> 
(1 ^ i ^ M ) 


It should be noted that the y. vectors do not form an orthog- 

0 

onal set. They can, however, be orthogonalized . That is, using 

— (2l 

the singular value decomposition for Y, V in (C.2) becomes 


V = EY = E(Ad^H') = GD^H' 


= GW = 




Ss^s 


(C.3) 


The rank of each of the matrices G, D^, and H is s where 

_ i 

s - r. Therefore 15^ has s positive elements along its diagonal 



and 


D® = G'VH (C.4) 


since H'H = G’G = 1^. 

Since the rank of the product of matrices is less than or 
equal to the rank of any of the individual factors, from (C.3) 
r ^ s and from (C.4) s - r. Therefore s = r. 

Letting S = W we see from (C.l) and (C.3) that there are 
two semi-orthogonal matrices E and G which diagonalize S. 


E'SE = YY' = D 


(C.5) 


G'SG = WW' = D 


(C.6) 


Since the characteristic equation is the same in both cases 
and the eigenvalues have been ordered (largest to smallest) D 
and D must be identical . Let D be the matrix . 


We must show now that the eigenvectors e. and g. are 

J J 

identical up to a sign change in their components. 


with 


Let e. and g. be two eigenvectors which correspond to d. 
J J J 


“ *®ij ®pj’ 


Sj = 


Then from (C.5) and (C.6) 






"*0 


which can be written as 



i,k 


i.k 


from which by equating like coefficients of the elements of S, 


X • 6 • f 


^ik’ 


we obtain 


'ij 


e, . e . . 
kj ij 


= g 


ij 


^kj%o 


(k =. i) 

(k + i) 


(C.7) 

(C.8) 


Equation (C.7) has two solutions 

= * Sij 

Substituting (C.9) into (C.8) gives 


®kj = i Skj 


and therefore the eigenvectors e- and g. are identical up to a 

J J 


sign change in their components. 


(C.9) 


The trace T of S can be written in terms of the trace of 


S- 

1 


1 1 


,M) 

r 

r = 

E 

0=1 


M 


i=l 


A \ 

. r: IE..,) ■ * 

i=l Vj=i 


. . .+ T, 


M 


(C.IO) 


(r^ + 


• • • 





where the (j 
the (J = 1,2 

elements of D- , 


= l,2,...,r) are the positive eigenvalues of S, 
,.,.,r^) are the positive eigenvalues of or 
and |a| is the Euclidean norm of A = (s-ij)* 


It should be noted that if i^,i 2 »...,ij 


is an arbitrary 

ordering of the integers from 1 to M, then there exists a (p x p) 
permutation matrix P such that PV will order the » • • • * 


the order V. 

^M, 


and the trace remains invariant. That is, 


trace (PV)(PV)' = trace (PV)’(PV) 


= trace (V’P'PV) = trace (V'V) = T 


(C.ll) 


Let us approximate V by selecting only eigenvectors (and 
corresponding principal components) from each subdivision such 
that lOOaCfi of the trace T is accounted for where 0 1. That 

is, let k^ - r^^ be the smallest integer such that 



(C.12) 


Then, using the first k^ terms in the singular decomposition 
of each V^, i.e., (see Appendix B) V is approximated by V 

(k = kj^ + . . . + kjyj) 


V ^ 


A 

V 






(C.13) 


The g and w vectors in (C.13) are the eigenvectors and components 
of VV . In (C.3) they are the eigenvectors and components of VV . If 
k^ “ “ 1,2,...,M), V = V and they then refer to the same matrix. 


Using (C.IO) and (C.12) 


(C.14) 


C 

n -1 



</T 


_ , A 

Now let k - k be 


the smallest integer such that 



- o^T 


where d • = w.w'. are the positive eigenvalues of YY' , 

J J J 

- dg - ... - d_ - d. > 0 and k^ ( i = 1,2, ... ,M) principal 
k k _ 

components from Y^ have been concatenated to form Y. 


The matrix V is then approximated by V 


^ o — 

V = V 


g^w^ + . . . + g_w_ 


k k 


(C.15) 


with 


A 

T 


= C<>i ^ Co. = 




j=l ^ 


f ^ o^T 


(C.l6) 


where T is the trace of YY' and T the trace of VV 


It should be noted that we have used V here as the best 

matrix in the least squares sense of rank k which approximates 

V, whereas in Appendix V was defined as the best matrix of 

rank k approximating V. If k. = r. for i = 1,2, ...,M, then V = V, 

A — ^ ^ 

k = r, k = k and the two definitions are consistent. In Appendix D 
for the first partitioning level (1=1) we set = V and then 

V is the least squares approximation to both V which is 

again consistent with Appendix B and this appendix. 



/I) 

' ' The partitioning and joining scheme described in this appendix 
was first presented by Yakovleva et. al (1968) and is documented 
in Rinne and Karhila (1979) • The difference between the presen- 
tation given here and that of Yakovleva is that the number of 
eigenvectors and principal components to be joined at the various 
levels is variable and depends upon the required tolerance level 
oC and the magnitudes of the largest eigenvalues within each 
subdivision. In the cited references a fixed number of vectors 
and components are joined together. 

The a. vectors (columns of A) or left sin^lar vectors of Y 
(alsoHhe eigenvectors of YY) are the "joining vectors" which 
are equivalent to the "joining functions" in Yakovleva. 


(2) 



APPENDIX D 


PARTITIONING AND JOINING IN GENERAL 


Let V "be a real (p x q) data matrix (q - p) with rank r 
where the row means have been subtracted from the data. Par- 
tition V into M subarrays with dimension (p^ x q) and rank 

r^ where p^^ + . . . + Pjyj = P. r^ + . . . + rjyj - r, r^ - q and 
r - q 



(D.l) 


From Appendix C if P is any (p x p) permutation matrix 
which reorders the 



(D.2) 


where i. (j = 1,2, are the integers i (i = 1,2 M) in 

J 

some order, then the trace remains unchanged. That is 


T = 





f 9 • 


+ 



(D.3) 


with T the trace of VV , 
Euclidean norm of A = 


^i 


the trace of 


(aij). 




Y. V! , and 


11 ' 

iTo 


|a| the 


Sets of local EOF's corresponding to the various subdivisions 
of the data array, V^, may be ordered and joined in a number of ways. 



The ordering of the M matrices in (D.2) is only 1 out of a possible 

Ml such permutations. For each permutation the V. can be joined to- 

M- 1 ^ M- 1 

gether in 2 possible combinations giving a total of M!2 

possible arrangements. In each case we are interested in retaining 

100 e( % of the total information content in the data (i.e., 100 o( % 

of the trace T) where is a selected tolerance level (0 ^ ^ 1). 


In general let there be 1 partitioning levels with matrices 
and groups of matrices at the 1 th level (1 = 1,2,...,L), and 
j^m^ matrices in the n th group which we will concatenate and join to- 


gether (Appendix C), 


l”l 


i “n 


= M 


1 


1 ’ 


The partitioning and joining procedure can now be described by 
the following matrices at the 1 th level 


h^(l) = ^(1) = \l+l) * 'f(l-M) 

(1 = 1,2, .. . ,L) 

where 



(D.4) 


(D.5) 


(D.6) 



(D.7) 


^( 1 + 1 ) 



The (elements of ) are least squares approximations 

A 

to the (i=l ,2, . . . ,M^) and consist of the first terms in the 

singular decomposition of (see Appendix C). is a permutation 

A _ A _ 

matrix. Initially, = M, = Ijj, “ ^(i) " ^i “ ^i ~ ^i 

for i = 1 , 2 s . . . sM. 


The elements of ere the permuted elements of which 
are grouped into groups, the nth group consisting of the foll- 
owing set of matrices 


n th group 






(D.8) 


(“n = “n-1 + “n- “o = 
(n 1,2,.. a, ^ 2 ^) 


where the 1 has been suppressed in ^m^ and ^m^ for simplicity. 


Joining the m^ matrices in (D.8) gives and element of 

A A 

V(i^l) (n = 1,2,...,N^). From Appendix C, can be written 

using the subscript n for joined matrices in the n th group 



(D.9) 


= E Y = E^(A D^H’) = 
n nn n nnn nnn nn 


The matrix consists of the joining vectors or functions 
which are also the left singular vectors of or equivalently 
the eigenvectors of Y^^Y^. 


Each of the m elements of E 
n n 



(D.IO) 


or the ^E^ (i = l,2,...,m^) matrices consists of a set of local 
EOF’s (columns of • 

In similar fashion each of the m elements of Y^ 

n n 


h\ 


(D.ll) 


or the ^Y^ principal component matrices consists of k^ principal 
components corresponding to the k^ EOF's in ^^E^. 

The columns of G^ are the "global" EOF's for the n th group 
formed from the ni sets of local EOF's. 



A _ 

Similar to Appendix C, for each the first k terms in 
its singular decomposition become the matrix in 

and we set k^ = ic and (n = 1,2 N^). 

From (D.4) 

= (Pl-.-Pi)V (D.12) 


and therefore 




(D.13) 


Denote by the trace of . We can write 


1 

W) 


(D.14) 


(1 = 1,2, . .. ,L) 


and then 

^(L+1) 


(D.15) 


We therefore can approximate V at a 100 of ^ level and each 
of the k^ is determined such that k^ = k where k is the 
smallest integer such that 



j=l 



^.1 



Table 1 


First 10 Eigenvalues 


Eigenvalue 

Number 

Eigenvalue 

(°K)-= 

% Contribution 
To Variance 

Cumulative 
Contribution ( 35 ^ 

1 

7026087.28 

68.78 

68.78 

2 

855772.20 

8.38 

77.16 

3 

458II7 .91 

4 .48 

81.64 

4 

320376.35 

3.14 

84.78 

5 

208065.12 

2.04 

86.82 

6 

142705.51 

1 .40 

88.21 

7 

111005.21 

1.09 

89.30 

8 

100980 .40 

0.99 

90.29 

9 

75789.78 

0.74 

91.03 

10 

61717.53 

0.60 

91.63 


Total variance = 10215181 (°K)^ 


> 



FIGURE 1: MAXIMUM ICE EXTENT AS REPRESENTED BY 
THE 155**K BRIGHTNESS ISOTHERM (FROM RAYNER AND 
HOWARTH) AND AREA ANALYZED 




FIGURE 2: MINIMUM ICE EXTENT AS REPRESENTED BY 
THE 155<*K BRIGHTNESS ISOTHERM (FROM RAYNER AND 
HOWARTH) AND AREA ANALYZED 




FIGURE 3: LOCATION MAP FOR SOUTH POLAR REGION 
(TAKEN FROM ZWALLY AND GLOERSEN 1977) SHOWING 

AREA ANALYZED 
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LEAST SQUARES POWER SPECTRUM 
OF THE SPATIAL MEAN 



LEAST SQUARn POWER SPECTRUM OF SPATIAL 
STANDARD DEVUTION 



FIGURE 4: SPATIAL MEAN AND STANDARD DEVIATION VERSUS TIME AND POWER SPECTRA 







Log of Eigenvalue 



FIGURE 5: LOGARITHM OF EIGENVALUE VERSUS 

NUMBER 



FIGURE 6: MEAN AND STANDARD DEVIATION MAPS, EOF’S 1 AND 2 AND CORRESPONDING PRINCIPAL 

COMPONENTS AND SPECTRA 
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FIGURE 7: EOF’S 3 AND 4, PRINCIPAL COMPONENTS AND SPECTRA 
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FIGURE 10: EOF’S 9 AND 10, PRINCIPAL COMPONENTS AND SPECTRA 







FIGURE 11: TIME HISTORY OF PIXEL 
BRIGHTNESS TEMPERATURE 







FIGURE 12: MONTHLY AVERAGE SEA ICE 
EXTENTS DERIVED FROM 1974 ESMR-5 IMAGES: 
(A) PERIOD OF ICE GROWTH; (B) PERIOD OF ICE 
DECAY (TAKEN FROM CAVALIERI AND 
PARKINSON) 




* Principal component r«lativt contribution it the percent contribution of the ith principal component to the total turn of squares of the jth measurement vector vj at the jth time point. 


FIGURE 13: PRINCIPAL COMPONENT RELATIVE CONTRIBUTION* 
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•PRINCIPAL COMPONENT RELATIVE CONTRIBUTION IS THE PERCENT CONTRIBUTION OF THE i"’ PRINCIPAL COMPONENT 
TO THE TOTAL SUM OF SQUARES OF THE j"* MEASUREMENT VECTOR y AT THE j*'' TIME POINT 


FIGURE 14: PRINCIPAL COMPONENT RELATIVE CONTRIBUTION* 






